function calculation_of_integral(CONSTS, plot_data)

    c     = CONSTS.c;
    k0    = CONSTS.k0;
    a     = CONSTS.a;
    d     = CONSTS.d;

    zeta_step = 0.05; % 0.05
    zeta_max = 1;
    zeta_vec = (zeta_step : zeta_step : zeta_max)'*d; 
    integ_vec = zeros(size(zeta_vec));
    integ_ch_min_vec = zeros(size(zeta_vec));
    integ_ch_max_vec = zeros(size(zeta_vec));
    for i = 1:size(zeta_vec, 1)
        zeta = zeta_vec(i);
        m = 1; % 0, 1, 2, 3, 4, ..., 20, ...
        period = pi/(k0*a);
        dots_per_period = 100; % 200
        order = 3;
        q_high = 0.5*log(10)*order/(k0*zeta);
        q_vec = (0.0 : period/dots_per_period : q_high)';

        y_vec_1 = func_1(q_vec, zeta, m, CONSTS); %zeros(size(q_vec));
        y_vec_2 = func_2(q_vec, zeta, m, CONSTS); %zeros(size(q_vec));
        integ_1 = trapz(q_vec, y_vec_1);
        integ_2 = trapz(q_vec, y_vec_2);

        integ_vec(i) = - ((2*pi*k0^2*a)/c) * integ_1 ...
                       + (2*pi/(c*a)) * integ_2;

        integ_ch_min_vec(i) = (2*k0/c) * integ_check_for_K1_1(zeta, m, CONSTS) + ...
                              (2/(c*k0*a^2)) * integ_check_for_K1_min_V(zeta, m, CONSTS);
        integ_ch_max_vec(i) = (2*k0/c) * integ_check_for_K1_1(zeta, m, CONSTS) + ...
                              (2/(c*k0*a^2)) * integ_check_for_K1_max_V(zeta, m, CONSTS);
    end   

    if (plot_data)
        figure; plot( ...
                     zeta_vec/d, integ_ch_min_vec, 'r+-', ...
                     zeta_vec/d, integ_ch_max_vec, 'rx-', ...
                     zeta_vec/d, integ_vec, 'b.-' ...
                     );
        title('K_m(\zeta)');
        xlabel('\zeta /d'); ylabel('K_m'); %ylim([0 8e-13]);
        legend('Symbolic integral min','Symbolic integral max', 'Numerical integral', 'Approximation', ...
            'location', 'SouthEast');
    end

end